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MRI  Assessment  of  the  Viscoelastic  Properties  of 
Normal  and  Abnormal  Breast  Tissue: 
Annual  Report  to  the  US  Army 


INTRODUCTION 

It  is  well  known  that  variations  in  tissue  elasticity  properties  are  associated  with 
the  presence  of  cancer.  Data  from  the  literature1,2  suggest  that  by  measuring  the  elastic 
properties  of  tumors  and  lesions  compared  to  normal  breast  tissue,  the  relationship 
between  pathological  processes  and  elasticity  can  be  investigated.  During  the  past 
several  years,  a  number  of  researchers3,4,5,6,7  have  been  attempting  to  use  MRI 
elastography  to  assess  the  biomechanical  properties  of  breast  tissue.  The  subject  of  this 
research  project  is  to  obtain  a  better  understanding  of  the  biomechanical  properties  of 
breast  tissue  by  accurately  measuring  their  stiffness  both  in-vivo  via  MRI  elastography 
and  invitro  using  a  precision  loading  system.  The  results  of  this  investigation  have  the 
potential  of  aiding  in  improved  diagnosis  of  breast  cancer  and  predicting  tissue 
deformation.  The  latter  is  very  important  in  various  medical  applications  such  image 
coregistration. 

Our  approach  to  measure  the  biomechanical  properties  of  breast  tissues  will  be 
to  apply  MRI  motion  detection  methods  to  measure  tissue  strain  while  the  surface  is 
undergoing  external  forces  with  an  applicator.  By  an  analysis  of  the  resulting  motions 
and  through  the  use  of  an  appropriate  elasticity  reconstruction  technique,  we  aim  to 
measure  the  distribution  of  tissue  elastic  properties  in-vitro  and  in-vivo.  These 
measurements  will  be  performed  on  a  number  of  normal  breast  tissue  samples  and  a 
number  of  breast  tissue  pathologies  will  be  obtained  from  pathology.  Based  on  these 
data,  the  feasibility  of  determining  the  type  of  breast  tissue  abnormality  and  predicting 
tissue  deformation  in-vivo  will  be  assessed. 


RESEARCH  TASK  -  Year  II 

This  project  is  a  three-year  project  with  a  number  of  separate  research  goals  and 
tasks.  During  the  second  year,  three  tasks  were  to  be  performed  as  follow: 

Task  I  -  The  first  task  is  to  test,  integrate  and  optimize  our  precision  table-top  loading 
system  for  in-vitro  elasticity  measurement  of  breast  tissues. 


Task  II  -  For  MRI  based  elasticity  measurement,  the  second  task  was  to  develop  a 
modulus  reconstruction  algorithm  which  inputs  the  MRI  derived  displacement  data  of  the 
sample  undergoing  quasi-static  loading  and  outputs  the  relative  elastic  moduli  of 
different  tissue  components  within  the  sample. 

Task  III  -  An  important  application  of  this  research  is  predicting  breast  tissue 
deformation.  Beside  the  fact  that  this  model  is  used  in  our  modulus  reconstruction 
algorithm,  it  can  be  applied  in  many  other  medical  applications  such  as  breast  image 
coregistration.  In  parallel  to  the  second  task,  the  third  task  was  to  develop  a  Finite 
Element  (FE)  model  that  requires  the  biomechanical  properties  of  breast  tissues  to 
predict  tissue  deformation.  This  model  has  been  successfully  used  as  the  forward  model 
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component  of  our  reconstruction  technique  as  well  as  to  demonstrate  other  potential 
medical  applications  of  this  project  such  as  breast  image  coregistration. 


PROGRESS  TO  DATE 

Task  I  -  Precision  table-top  system 

During  the  first  year,  we  started  the  design  and  configuration  of  a  table-top  uniaxial 
system  to  measure  tissue  elastic  properties  in-vitro.  Throughout  the  second  year,  this 
system  has  been  tested,  modified  and  integrated  to  provide  both  accurate 
measurements  as  well  as  conditions  similar  to  those  of  in-vivo.  This  system  is  shown  in 
Figure  1  and  it  allows  for  table-top  uniaxial  loading  experiments  under  a  broad  range  of 
frequencies  and  amplitudes.  This  system  consists  of  two  main  components:  a  linear 
stepper  motor  actuator  used  to  deform  the  tissue  sample  in  a  programmable  fashion  and 
a  load  cell  for  force  measurement.  The  applied  tissue  deformation  is  measured  based  on 
the  motion  delivered  by  the  actuator,  while  measuring  the  resulting  force  dynamically 
with  the  load  cell.  The  actuator  and  load  cell  are  wired  to  a  computer  that  gives 
instructions  to  the  actuator  and  receives  data  from  both  the  actuator  and  the  load  cell 
through  Labview  virtual  instrument  software.  To  provide  conditions  similar  to  those  of  in- 
vivo,  tissue  samples  are  kept  moist  by  immersing  them  in  an  isotonic  sodium  chloride 
solution,  and  temperature  is  maintained  at  37°c  using  a  temperature  control  device.  This 
temperature  control  device  is  yet  to  be  tested  with  biological  tissues.  The  loading  system 
has  been  tested  with  chicken  breast  tissue.  Five  samples  were  prepared  and  their 
Young’s  modulus  measured.  As  a  result,  the  mean  value  and  standard  deviation  of 
these  measurements  were  5.3  Kpa  and  0.4  kPa  respectively  which  suggests  that  the 
system  provides  both  accurate  and  repeatable  measurements. 

Task  II  -  Constrained  Modulus  Reconstruction  Technique 

During  the  first  year  of  this  project,  an  MR  compatible  loading  system  was  designed  and 
configured  which  allows  controlled  compression  of  small  tissue  samples  and  phantoms 
for  testing  purposes.  In  addition,  Navier  equations  based  moduli  reconstruction 
algorithms  were  developed  to  be  used  in  conjunction  with  the  loading  system.  To 
achieve  a  reasonable  reconstruction  results,  however,  very  high  SNR  values  are 
required  which  are  very  difficult  to  obtain  in  practice.  As  an  alternative,  we  have 
developed  an  iterative  constrained  modulus  reconstruction  technique  with  moderate 
SNR  requirement.  This  technique  has  been  tested  on  a  multi-component  phantom  with  a 
3-D  complex  geometry  and  the  results  indicate  that  the  reconstruction  algorithm  is 
reasonably  accurate  and  robust.  Figure  2  depicts  the  phantom’s  geometry  and  the 
modulus  reconstruction  results  obtained  from  this  technique.  A  manuscript  outlining  the 
mathematics  and  results  of  this  technique  has  been  submitted  for  publication  and  is  now 
in  press  in  lEEE-Transactions  on  Medical  Imaging.  To  apply  this  technique  for  in-vivo 
breast  tissue  modulus  reconstruction,  a  breast  compression  system  shown  in  Figure  3 
has  been  designed  and  constructed.  This  system  consists  of  two  compression  plates, 
one  stationary  and  the  other  is  driven  by  the  MR  compatible  stepper  motor  to  compress 
the  breast  sinusoidally.  At  this  stage,  this  system  is  being  tested  on  a  breast  shaped 
phantom  shown  in  Figure  4.  This  figure  shows  a  sagittal  MR  image  of  this  phantom  with 
a  preliminary  coronal  strain  image  of  the  phantom.  Although  encouraging,  the  strain 
image  suggests  that  further  improvement  to  the  compression  system  is  required  to 
achieve  better  SNR. 
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Task  III  -  Breast  Tissue  Deformation  Model 


Modeling  breast  tissue  deformation  is  a  very  important  task  for  this  research  project  as  it 
is  the  primary  component  of  our  constrained  elasticity  reconstruction  technique.  For  this 
reason,  in  parallel  to  the  modulus  reconstruction  technique,  we  have  developed  3-D  FE 
model  that  uses  breast  MRI  images  and  tissue  mechanical  properties  as  an  input  and 
calculates  tissue  deformation  and  stress  distribution  as  an  output.  This  model  is  capable 
of  predicting  both  small  deformation  as  in  the  case  of  breast  elastography  as  well  as 
large  deformation  as  required  in  other  applications  such  as  image  coregistration.  This, 
however,  requires  the  hyperelastic  properties  of  the  breast  tissues  to  be  known.  The 
details  of  this  model  and  its  application  in  breast  image  coregistrationis  presented  in  a 
paper8  which  was  published  recently  in  lEEE-Transactions  on  Medical  Imaging. 


FUTURE  WORK 

The  research  that  was  originally  proposed  for  this  project  remains  valid.  Conclusions 
from  the  first  and  second  years’  research  suggest  that  it  is  feasible  to  utilize  MRI  imaging 
to  measure  breast  tissue  biomechanical  properties.  This  concept  has  been  proven  using 
experimental  data  with  a  phantom  and  simulated  data  with  a  volunteer’s  breast  as 
presented  in  our  modulus  reconstruction  paper8.  After  testing  the  temperature  control 
system  mounted  on  our  table-top  uniaxial  loading  system,  we  will  start  measurements  of 
breast  tissue  elastic  modulus  in-vitro.  In  the  mean  time,  after  improving  our  breast 
compression,  we  will  conduct  in-vivo  breast  elastography  experiments  in  conjunction 
with  our  modulus  reconstruction  algorithm.  In  addition,  we  will  attempt  to  measure  breast 
hyperelastic  properties  in-vitro  which  are  very  important  in  predicting  breast  tissue  large 
deformation. 


KEY  RESEARCH  ACCOMPLISHMENTS 

1)  The  construction  of  a  precision  table-top  loading  system  to  perform  in-vitro  breast 
tissue  elasticity  measurement.  This  system  has  been  tested  and  proven  to  be 
relatively  accurate  and  reliable.  In-vitro  breast  tissue  elasticity  measurement  will  be 
started  in  the  near  future. 

2)  The  development  of  a  constrained  modulus  reconstruction  technique  for  in-vivo  breast 
elastography  using  MRI  derived  displacement  data.  A  paper  outlining  this  technique 
will  be  published  shortly  in  lEEE-Transactions  on  Medical  Imaging.  This  paper 
suggests  that  very  moderate  SNR  values  are  required  for  reasonably  accurate 
modulus  reconstruction.  In  addition,  a  breast  compression  device  has  been 
constructed  and  tested  with  a  breast  shaped  phantom.  The  preliminary  results  of  this 
test  is  encouraging,  however,  further  improvement  is  required  to  achieve  better  SNR 
values  required  for  in-vivo  experiments. 

3)  The  development  of  a  breast  tissue  deformation  model  using  MRI  data  has  been 
completed  and  a  paper  outlining  the  mathematics  and  method  has  been  published 
recently  in  lEEE-Transactions  on  Medical  Imaging.  This  model  is  capable  of 
predicting  the  displacements  and  stress  distribution  of  breast  tissues  while 
undergoing  both  small  and  large  deformation. 
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REPORTABLE  OUTCOMES 


During  the  past  year,  we  have  submitted  a  number  of  papers  and  symposium  abstracts, 
as  detailed  below,  which  relate  to  the  research  contained  in  this  project. 

PUBLISHED  OR  SUBMITTED  MANUSCRIPTS: 

1)  Samani  A,  Bishop  J,  Yaffe  MJ  and  Plewes  DB.  Biomechanical  3-D  finite  Element 
Modeling  of  the  Human  Breast  Using  MRI  Data.  IEEE  Transactions  on  Medical 
Imaging  20(4):271-279,  2001. 

2)  Samani  A,  Bishop  J  and  Plewes  DB.  A  Constrained  Modulus  Reconstruction 
Technique  for  Breast  Cancer  Assessmant.  IEEE  Transactions  on  Medical 
Imaging,  2001  (in  press). 

3)  Bishop  J,  Samani  A,  Sciarretta  J,  Luginbhul  C  and  Plewes  DB.  A  Signal/Noise 
Analysis  of  Quasi-static  MR  Elastography.  Submitted  to  IEEE  Transactions  on 
Medical  Imaging,  2001 

ABSTRACTS: 

1 )  Samani  A,  Bishop  J  and  Plewes  D.  A  Constrained  Breast  Magnetic  Resonance 
Elastography  Technique:  3D  Phantom  Study.  International  Society  for  Magnetic 
Resonance  in  Medicine:  Nineth  Annual  Meeting,  Glasgow,  April  21-27,  2001, 
Paper#1640. 

2)  Samani  A,  Bishop  J  and  Plewes  D.  3D  Finite  Element  Model  for  Breast  Non-rigid 
Registration.  International  Society  for  Magnetic  Resonance  in  Medicine:  Nineth 
Annual  Meeting,  Glasgow,  April  21-27,  2001,  Paper#837. 

3)  Bishop  A,  Samani  A  and  Plewes  D.  A  Signal/Noise  Analysis  of  Magnetic 
Resonance  Strain  Imaging.  International  Society  for  Magnetic  Resonance  in 
Medicine:  Nineth  Annual  Meeting,  Glasgow,  April  21-27,  2001,  Paper#1646. 


CONCLUSIONS 

During  the  first  and  second  years  of  this  project,  we  have  made  good  progress 
toward  our  end  goal  of  understanding  the  biomechanical  properties  of  breast  tissues  and 
developing  the  means  to  apply  these  in  tissue  deformation  analyses.  The  overall 
objective  of  this  project  was  to  create  reliable  mathematical  tools  to  predict  tissue 
deformation  in  the  presence  of  compression  to  determine  tissue  properties  of  various 
breast  tissues,  both  normal  as  well  as  benign  and  neoplastic  tissues.  We  have,  thus  far, 
been  able  to  develop  the  hardware  and  software  tools  required  to  conduct  this  study.  In 
addition,  we  have  used  these  tools  with  phantoms  and  animal  tissues  and  the  results 
indicate  that  they  are  reasonably  accurate  and  reliable.  During  the  third  year  of  this 
project,  we  will  use  these  tools  to  measure  the  biomechanical  properties  of  breast 
tissues  using  samples  we  will  obtain  from  pathology.  The  importance  of  this  work 
remains  high  providing  a  platform  for  new  detection  and  diagnosis  methods  for  breast 
cancer.  This  work  can  also  open  up  new  applications  of  breast  imaging  including  data 
fusion  and  therapy  planning  tools  based  on  predictive  tissue  deformation  analysis 
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capabilities  such  as  the  breast  tissue  deformation  model  we  developed  in  the  second 
year  of  this  project. 
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Appendix  1 :  Figures  of  the  body  of  the  proposal 


Figure  1.  Precision  table-top  uniaxial  loading  system  for  in-vitro  tissue  measurement.  This 
system  consists  of  two  main  component:  a  linear  stepper  motor  actuator  used  to  deform  the 
tissue  sample  in  a  programmable  fashion  and  a  load  cell  for  force  measurement. 


Figure  2.  MRI  image  of  the  multi-layer  phantom  (left)  and  reconstruction  results 
using  the  constrained  reconstruction  technique  (right).  The  reconstruction  results 
indicate  a  maximum  error  of  12%  in  reconstructing  the  modulus  of  different  layers 
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Figure  3.  MRI  compatible  breast  compression  system  used  for  MR  elastography 
of  a  breast-shaped  phantom.  This  system  consists  of  commercial  breast  coil  and 
an  assembly  of  plates  and  cams  which  deliver  sinusoidal  compression  to  the 
breast  via  an  MR  compatible  stepper  motor. 


Figure  4.  MRI  sagittal  image  of  the  breast-shaped  phantom  (left)  and  a  strain 
image  a  coronal  slice  corresponding  to  section  a-a  (right).  Although  the  strain 
SNR  is  not  satisfactory,  different  layers  of  the  phantom  are  visible. 
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Appendix  2:  Published  abstracts  and  papers 
Abstracts: 

1)  Samani  A,  Bishop  J  and  Plewes  D.  A  Constrained  Breast  Magnetic  Resonance 
Elastography  Technique:  3D  Phantom  Study.  International  Society  for  Magnetic 
Resonance  in  Medicine:  Nineth  Annual  Meeting,  Glasgow,  April  21-27,  2001, 
Paper#1640. 

2)  Samani  A,  Bishop  J  and  Plewes  D.  3D  Finite  Element  Model  for  Breast  Non-rigid 
Registration.  International  Society  for  Magnetic  Resonance  in  Medicine:  Nineth 
Annual  Meeting,  Glasgow,  April  21-27,  2001,  Paper#837. 

3)  Bishop  A,  Samani  A  and  Plewes  D.  A  Signal/Noise  Analysis  of  Magnetic 
Resonance  Strain  Imaging.  International  Society  for  Magnetic  Resonance  in 
Medicine:  Nineth  Annual  Meeting,  Glasgow,  April  21-27,  2001,  Paper#1646. 
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k  Constrained  Breast  Magnetic  Resonance  Elastography  Technique:  3D  Phantom  Study 


A.  Samani,  J.  Bishop,  and  D.  B.  Plewes 
Department  of  Medical  Biophysics,  University  of  Toronto 


Introduction  Tissue  elasticity,  defined  by  Young’s  modu¬ 
lus,  is  well  known  to  be  associated  with  the  presence  of  cancer. 
Compared  to  normal  fibroglandular  tissue,  recent  studies  [1] 
suggested  a  6  fold  increase  in  the  Young’s  modulus  of  breast 
carcenoma  while  fibroadenomas  are  only  a  factor  of  3  stiffer. 
These  observations  indicate  that  an  imaging  modality  which  is 
based  on  elasticity  contrast  can  improve  specificity.  This  have 
stimulated  the  development  of  breast  MR  elastography  (MRE) 
techniques  [2,3]  with  the  aim  of  imaging  tissue  Young’s  modu¬ 
lus.  In  these  techniques,  harmonic  [2]  or  quasi-static  [3]  motion 
is  mechanically  induced  within  the  tissue,  the  resulting  tissue 
deformation  is  measured  by  MRI  techniques,  from  which  the 
tissues’  Young’s  modulus  distribution  can  be  calculated  via  an 
inversion  technique.  In  [3],  we  presented  a  quasi-static  MRE 
technique  which  proved  to  be  successful  based  on  a  simple  2d 
phantom  study,  and  a  breast  study  using  simulated  data.  In 
this  article,  we  present  the  results  of  a  3d  phantom  consisting 
of  three  different  layers  with  complex  geometry  to  simulate  adi¬ 
pose,  fibroglandular  and  tumour  tissues.  Also,  for  breast  MRE, 
modulus  reconkruction  results  using  an  efficient  finite  element 
(FE)  technique  is  presented. 

Methods  This  technique  is  similar  to  the  iterative  recon¬ 
struction  MRE  technique  presented  in  [3]  where  the  geometry 
of  each  tissue  type  is  assumed  to  be  known  apriori .  Quasi¬ 
static  sinusoidal  compression  of  0-5  mm  amplitude  is  applied 
at  the  surface  of  the  phantom  using  a  MR  compatible  device 
and  the  resulting  deformation  is  measured  using  a  stimulated 
echo  (STEAM)  pulse  sequence.  The  iterative  modulus  recon¬ 
struction  process  consists  of  tissue  stress  calculation  using  the 
.  FE  method  followed  by  modulus  updating  of  each  element  using 
Hooke’s  law.  At  each  iteration,  the  modulus  of  each  tissue  type 
is  calculated  by  averaging  over  its  corresponding  volume.  For 
breast  MRE,  we  used  a  transfinite  interpolation  (TI)  technique 
[4]  which  is  based  on  the  idea  of  mapping  a  unit  square  to  any 
given  enclosed  region.  This  technique  can  model  the  breast’s 
curved  surface  more  efficiently  than  the  previously  used  voxel 
based  meshing  technique.  Moreover,  this  technique  is  capable 
of  using  variable  FE  size  which  makes  it  possible  to  increase 
the  number  of  elements  involved  in  averaging  in  the  tumour,  if 
necessary. 

Results  A  cubic  plastisol  PVC  phantom  composed  of  a 
a  12  mm  spherical  stiff  volume  representing  a  tumour  sur¬ 
rounded  by  a  softer  volume  representing  fibroglandular  tissue 
was  constructed.  The  latter  was  enclosed  by  another  softer 
third  layer  representing  adipose  tissue.  The  Young’s  moduli 
of  the  individual  layers  were  measured  independently  using 
a  uniaxial  benchtop  experiment,  resulting  in  48.0  kPa,  24.3 
kPa  and  12.4  kPa  for  the  three  layers.  MRI  images  of  this 
phantom  was  acquired  using  a  3d  fast  spin-echo  sequence  with 
TR/TE=2000/85  ms  and  48  echos.  The  phantom  was  placed 
in  the  compression  device  where  it  underwent  a  quasi-static 
sinusoidal  compression  of  3.5  mm  amplitude.  A  STEAM  pulse 
sequence  with  TR/TE=660/12  ms  and  a  mixing  time  Tm  of 
300  ms  was  used  to  acquire  displacement  data  in  the  central 
plane  shown  in  Figure  1.  Starting  with  Young’s  modulus  ratios 
obtained  from  the  measured  strain  image  as  an  initial  guess, 
as  shown  in  Figure  1,  convergence  was  achieved  after  11  iter¬ 
ations.  The  final  reconstructed  modulus  ratios  indicate  that 
there  is  a  maximum  error  of  12%.  For  breast  MRE,  we  used 


a  volunteer’s  breast  MRI  image  to  create  a  3d  FE  mesh  using 
the  TI  based  method.  A  sagittal  slice  of  this  mesh,  where  a 
simulated  tumour  is  added,  is  depicted  in  Figure  2. 


Figure  1:  a)  FE  mesh  of  the  phantom’s  central  plane  and  b) 
reconstructed  modulus  ratios  where  the  dashed  lines  represent 
the  measured  modulus  ratios. 

The  Young’s  modulus  of  the  adipose,  fibroglanular,  tumour  and 
skin  tissues  were  assumed  to  be  2.0  kPa,  10.0  kPa,  50.0  kPa 
and  25.0  kPa  respectively.  Accordingly,  displacements  resulting 
from  5  mm  compression  normal  to  the  sagittal  plane  were  cal¬ 
culated  using  a  contact  problem  FE  model.  The  displacement 
component  in  the  compression  direction  was  contaminated  by 
normally  distributed  noise  to  simulate  SNR  value  of  10.  Using 
this  displacement  component  and  starting  with  modulus  ratios 
obtained  from  the  strain  image,  as  shown  in  Figure  2,  conver¬ 
gence  to  the  true  values  was  achieved  after  7  iterations. 
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Figure  2:  a)  FE  mesh  of  a  breast  slice  b)  MRE  results  where 
the  dashed  fines  represent  the  exact  modulus  ratios. 

Conclusions  This  article  presents  new  developments  in 
our  previous  constrained  breast  MRE  technique  [3].  This  tech¬ 
nique  is  very  efficient  as  it  requires  imaging  only  one  displace¬ 
ment  component  in  one  slice.  Using  a  phantom  consisting  of  3 
layers  with  complex  3D  geometry,  it  was  shown  that  this  tech¬ 
nique  is  capable  of  solving  the  breast  MRE  problem  where  the 
3D  stress  and  strain  calculation  as  well  as  image  segmentation 
are  challenging.  Furthermore,  to  improve  the  performance  of 
our  technique,  we  used  a  TI  based  FE  meshing  technique  which 
leads  to  better  geometry  representation  and  allows  using  vari¬ 
able  element  size  in  the  FE  mesh.  The  latter,  as  demonstrated 
in  the  breast  MRE  example,  can  lead  to  a  larger  number  of 
elements  in  the  tumour  to  suppress  the  noise  effects  through  av¬ 
eraging  while  the  FE  problem  size  does  not  change  significantly. 
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Introduction  Various  medical  applications  often  require 
comparison  or  data  fusion  of  MR  images  that  are  acquired  at 
different  times.  For  example,  to  determine  a  breast  tumour 
growth  over  time,  radiologists  need  to  compare  the  latest  MR 
image  with  MR  images  acquired  previously.  However,  since  it  is 
inevitable  that  MR  images  are  acquired  under  different  breast 
compressions  or  body  positions,  tissues  will  deform  and,  thus, 
lead  to  significant  changes  in  the  anatomical  geometry.  As  such, 
image  registration  is  required  to  facilitate  image  comparison  or 
data  fusion.  Many  researchers  such  as  [1]  have  used  voxel-based 
similarity  measure  concepts  in  their  registration  algorithms  ig¬ 
noring  the  physics  behind  tissue  deformation.  Based  on  an 
alternate  methodology,  in  this  article,  we  present  a  non-rigid 
image  registration  model  for  the  breast  based  on  biomechanics. 
We  use  a  3D  nonlinear  Finite  Element  (FE)  model  as  a  numer¬ 
ical  solution  to  solve  the  breast  tissue  deformation  problem. 

Methods  and  Theory  We  use  a  3D  FE  model  based  on 
the  finite  deformation  theory  of  elasticity  in  conjunction  with 
hyperelastic  material  to  model  large  deformations  of  the  breast 
undergoing  medical  procedures.  In  a  static  condition,  the  equi¬ 
librium  equation  of  an  object  under  applied  forces  is: 


Figure  2a,  were  acquired  with  a  GE  SIGN  A  1.5  T  scanner.  We 
used  a  2-D  spin-echo  pulse  sequence  with  TR/TE  =  300/9  ms 
and  9  =  90°.  The  objective  here  was  to  calculate  MR  images 
of  this  breast  as  if  it  underwent  8.0  mm  lateral  compression 
by  two  rigid  plates.  For  this  purpose,  a  FE  mesh  was  created 
and  the  compression  was  simulated  using  a  3D  contact  problem 
model.  A  FE  mesh  of  a  typical  slice  and  a  schematic  of  the 
contact  problem  model  are  depicted  in  Figure  1.  For  the  fat 
and  fibroglandular  tissues,  Neo-Hookean  hyperelastic  models 
were  calculated  using  Wellman  1999  data  [4]  while  the  skin  was 
assumed  to  be  elastic  with  a  Young’s  modulus  of  25.0  kPa. 


a) 

Figure  1:  a)  FE  mesh  of  a  sagittal  slice  and  b)  contact  prob¬ 
lem  model  of  the  breast. 


9  (Tn  dtTj  2 
dx-L  8X2 


9(Tj3 

dx3 


+  /i  =  0; 


i  — 1,2,3  (1) 


where  the  tr’s  denote  the  Cauchy  (true)  stress  tensor  compo¬ 
nents  and  fi  represents  the  body  forces.  Assuming  that  a  strain 
energy  function  U  is  selected  as  a  function  of  the  strain  invari¬ 
ants  Ji,  I2  and  J,  a  can  be  obtained  from  [2]: 


♦s>  « 


where  B  is  a  function  of  the  deformation  gradient  and  I  is  the 
identity  matrix.  We  use  ABAQUS:  a  commercial  FE  software 

[2]  which  solves  equations  (1)  and  (2)  numerically  to  calculate 
the  displacements.  For  FE  mesh  generation  using  hexahedral 
elements,  MR  images  are  first  segmented  using  AnalyzeAVW 
2.5  to  detect  breast  edges  and  separate  fat  and  fibroglandu¬ 
lar  tissues .  The  segmented  images  are,  thus,  processed  using 
a  transfinite  interpolation  meshing  technique  as  described  in 

[3] .  Assuming  that  the  displacement  boundary  conditions  are 
known,  the  breast  tissue  displacement  field  can  be  calculated. 
This  displacement  field  in  conjunction  with  the  MR  image  of 
the  breast  acquired  prior  to  the  movement  is  used  to  calculate 
an  MR  image  consistent  with  the  new  configuration. 


Results  The  mesh  generation  technique  was  first  validated 
using  a  60  x 60  x  60  mm  cubic  phantom  made  of  agar  with  a  hard 
cylindrical  inclusion.  MR  image  of  the  phantom  was  acquired 
and  a  FE  mesh  was  created  after  image  segmentation.  Another 
FE  mesh  was  created  manually  based  on  the  phantom’s  known 
geometry.  Using  these  meshes,  a  lateral  compression  of  9.0 
mm  was  simulated  with  ABAQUS  and  the  corresponding  dis¬ 
placements  were  calculated  at  the  inclusion’s  interface  where 
maximum  error  is  expected.  As  a  result,  the  average  difference 
between  these  displacements  was  less  than  0.5%  For  breast 
image  nonrigid  registration,  Ti  weighted  sagittal  MR  images 
of  a  healthy  volunteer,  a  typical  image  of  which  is  shown  in 


Using  ABAQUS,  the  problem  was  solved  and  the  displacements 
were  calculated.  These  displacements,  in  conjunction  with  the 
precompresaion  MR  image,  were  used  to  calculate  the  post 
compression  MR  image.  A  sagittal  image  and  its  correspond¬ 
ing  calculated  post  compression  image  are  shown  in  Figure  2. 


Figure  2:  a)  Precompression  MR  image  and  b)  calculated 
post  compression  image  of  the  breast. 


Conclusions  This  article  presents  a  biomechanics  based 
3D  FE  model  which  can  be  used  for  breast  image  nonrigid  reg¬ 
istration.  This  model  is  based  on  the  finite  deformation  theory 
of  elasticity.  For  mesh  generation,  we  used  a  meshing  technique 
we  have  developed  using  the  transfinite  interpolation  method. 
This  technique  was  validated  and  found  to  be  accurate  for 
displacement  calculation.  Breast  tissues  were  modeled  using 
hyperelastic  Neo-Hookean  models  while  the  skin  was  assumed 
to  be  elastic.  "The  results  shown  in  Figure  2,  qualitatively, 
prove  the  merits  of  this  model.  However,  further  research  is 
required  to  test  the  model  quantitatively. 
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Introduction  In  MR  elastography,  images  of  displacement  within 
tissue  or  a  phantom  are  obtained  during  a  repeatable  deformation. 
We  and  others  have  used  stimulated-echo  phase-contrast  MR 
imaging1,2,3  to  obtain  the  displacement  data.  Strain,  which  is  the 
spatial  gradient  of  displacement,  represents  a  simple  but  effective 
means  of  displaying  die  results  of  MR  elastography  because  it  is 
easy  to  compute  and  relatively  straightforward  to  interpret  In  this 
abstract,  we  extend  the  analysis  of  Steele  et  al4  by  describing  and 
validating  an  analytic  expression  for  strain/noise  ratio.  The 
expression  is  applicable  to  stimulated  echo  methods  and  can  be 
used  to  design  an  elastography  experiment  for  maximum  strain 
SNR. 


Theory  Assume  that  a  uniform  object  is  compressed  from  the 
top  surface  in  the  vertical  (y)  direction  and  is  fixed  at  the  bottom 
surface.  If  the  pulse  sequence  of  figure  1  is  used  to  measure  the  y 
displacement  by  phase-contrast,  then  the  phase  signal  within  the 
object  ©  will  be  proportional  to  vertical  displacement  u  as  a 
function  of  height  y: 

Q(y)=2n&du(y),  ©(0)=k(0)=  0  (1) 

where  d>d=yGx  is  the  zeroth  moment  of  the  displacement  encoding 
gradient  pulse.  If  L  is  the  height  of  the  object,  then  the  strain  is 
given  by  the  average  phase  gradient  is: 

30  =  ©(L)_£ 

dy  ~  L  h  U 

where  0  is  the  average  phase  dispersion  per  pixel,  and  h  is  the  size 
of  a  pixel. 


The  well-known  relation  for  phase  noise  per  pixel  is 


where  SNR  is  the  signal/noise  of  the  magnitude  data.  In  computing 
the  discrete  phase  gradient  of  (2)  by  a  central  difference 
approximation,  the  phase  noise  becomes: 


1  1 
-J2h  -J2SNR 


(4  . 


where  h  is  the  pixel  size  from  (2).  The  ratio  of  (2)  and  (4)  can  then 
be  interpreted  directly  as  the  strain/noise  ratio: 

SNRslrain  =2SNR8,  Q  <n  (5) 


If  0  >  tc,  the  imaging  point-spread  function  eliminates  all  signal5. 
Finally,  the  magnitude  SNR  term  in  (5)  is  expanded  by  explicitly 
writing  the  magnitude  signal  attenuation  caused  by  the  diffusion 
sensitivity  of  the  displacement  encoding  pulses: 

SNRstrain=2SNRe  3  0,  Q  <n  (6) 

where  D  is  the  apparent  diffusion  coefficient. 

Methods  Imaging  experiments  were  conducted  with  quasi-static 
compression  and  a  stimulated  echo  pulse  sequence  (Figure  1)  gated 
to  the  compression  cycle.  A  uniform  cubic  phantom  of  dimension 
25  mm,  made  of  1 %  agar,  was  compressed  by  1  mm  at  a  frequency 
of  2  Hz..  Imaging  time  was  adjusted  to  achieve  a  base  SNR  (ie  with 
no  displacement  encoding)  in  the  magnitude  data  of  about  30.  Data 
were  collected  for  a  range  of  displacement-encoding  sensitivities  at 
constant  echo  time.  Mean  strain  and  strain  noise  were  measured 
from  experimental  data  and  compared  to  the  predictions  of  (6). 


Results  Figure  2  plots  the  theory  ( - )  and  experimental 

results  for  both  the  magnitude  (o)  and  strain  (+)  SNR.  The 
encoding  sensitivity  is  plotted  on  the  x-axis,  expressed  in  terms  of 
the  phase  dispersion  per  pixel  0.  Note  that  0  and  encoding 
sensitivity  <t>d  are  related  through  the  surface  displacement  u(L)  in 
(1)  and  (2).  The  SNR  of  the  magnitude  data  drops  off  sharply  as  0 
increases,  due  to  diffusion  sensitivity.  Although  the  diffusion 
coefficient  was  not  independently  measured,  a  value  of 
D=2.05xl0'9  m2/s  fits  the  data  well. 

As  expected  the  strain  SNR  (+)  initially  increases  with  increasing 
amount  of  phase  dispersion  0.  The  strain  SNR  goes  through  a 
maximum  value,  then  quickly  drops  off  to  zero  well  before  the 
theoretical  maximum  encoding  of  n  radians  per  pixel  can  be 
achieved,  indicating  that  the  experiment  is  limited  by  diffusion 
sensitivity.  The  predicted  strain  SNR  from  (6)  is  shown  by  the 
solid  line,  with  good  agreement  to  experiment 

Conclusions  We  have  derived  and  validated  a  simple 
expression  for  strain  SNR  in  quasi-static  MR  elastography  using 
stimulated  echoes.  The  presentation  separates  the  diffusion 
sensitivity  from  all  other  terms  contributing  to  magnitude  SNR, 
and  simplifies  the  experimental  process  of  choosing  the  encoding 
sensitivity.  For  small  surface  displacements,  the  diffusion 
sensitivity  of  the  encoding  pulses  clearly  limits  the  achievable 
strain  SNR.  Although  stimulated  echoes  are  preferable  for  static 
and  quasi-static  methods  where  velocities  are  small,  spin-echoes 
can  be  used  to  reduce  diffusion  sensitivity  at  the  possible  expense 
of  increased  T2  sensitivity  in  longer  echo  times. 
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.  Biomechanical  3-D  Finite  Element  Modeling  of  the 

Human  Breast  Using  MRI  Data 

Abbas  Samani*,  Jonathan  Bishop,  Martin  J.  Yaffe,  and  Donald  B.  Plewes 


Abstract — Breast  tissue  deformation  modeling  has  recently 
gained  considerable  interest  in  various  medical  applications.  A 
biomechanical  model  of  the  breast  is  presented  using  a  finite 
element  (FE)  formulation.  Emphasis  is  given  to  the  modeling  of 
breast  tissue  deformation  which  takes  place  in  breast  imaging 
procedures.  The  first  step  in  implementing  the  FE  modeling 
(FEM)  procedure  is  mesh  generation.  For  objects  with  irregular 
and  complex  geometries  such  as  the  breast,  this  step  is  one  of 
the  most  difficult  and  tedious  tasks.  For  FE  mesh  generation, 
two  automated  methods  are  presented  which  process  MRI  breast 
images  to  create  a  patient-specific  mesh.  The  main  components 
of  the  breast  are  adipose,  fibroglandular  and  skin  tissues.  For 
modeling  the  adipose  and  fibroglandular  tissues,  we  used  eight 
noded  hexahedral  elements  with  hyperelastic  properties,  while  for 
the  skin,  we  chose  four  noded  hyperelastic  membrane  elements. 
For  model  validation,  an  MR  image  of  an  agarose  phantom  was 
acquired  and  corresponding  FE  meshes  were  created.  Based 
on  assigned  elasticity  parameters,  a  numerical  experiment  was 
performed  using  the  FE  meshes,  and  good  results  were  obtained. 
The  model  was  also  applied  to  a  breast  image  registration  problem 
of  a  volunteer’s  breast.  Although  qualitatively  reasonable,  further 
work  is  required  to  validate  the  results  quantitatively. 

Index  Terms — Breast,  finite  element  analysis,  image  registration, 
mesh  generation,  MRI. 


I.  Introduction 

PREDICTING  breast  tissue  deformation  is  of  great  sig¬ 
nificance  in  several  medical  applications  such  as  surgery, 
biopsy  and  imaging.  In  breast  surgery,  surgeons  are  often  con¬ 
cerned  with  a  specific  portion  of  the  breast,  e.g.,  tumor,  which 
must  be  located  accurately  beforehand.  Locating  the  portion 
of  interest  is  done  using  imaging  techniques  which  are  done 
under  tissue  configuration  of  compression  or  body  position 
that  is  often  entirely  different  from  that  of  surgery.  In  this  case, 
nonrigid  image  registration  is  required  to  locate  the  portion  of 
interest  in  the  new  configuration.  In  breast  biopsy,  the  biopsy 
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needle  causes  the  breast  tissue  to  deform  which  may  lead 
the  tumor  to  displace.  This  deformation  renders  the  guiding 
images  less  accurate  unless  the  displacements  are  determined 
and  needle  aiming  is  corrected  accordingly.  In  multimodality 
imaging  of  the  breast,  images  of  different  modalities  are  often 
obtained  under  different  tissue  configurations  of  compression, 
orientation,  or  body  position.  In  this  case,  determining  the 
tissue  deformation  is  essential  for  facilitating  data  fusion  [1], 
FEM  has  been  widely  used  in  modeling  biological  tissues  such 
as  bone  [2]~[4]  and  myocardium  [5],  and  recently  in  modeling 
brain  deformation  [6]-[8].  This  work  is  aimed  at  using  FEM  in 
modeling  breast  tissue  deformation. 

For  a  FE  model  to  be  satisfactory  in  predicting  deformation 
and  other  relevant  parameters,  in  addition  to  using  accurate 
elasticity  parameters  and  boundary  conditions,  a  reasonably 
accurate  FE  geometry  model  is  essential.  As  the  breast  is 
composed  of  soft  tissues,  large  deformation  is  usually  created 
while  undergoing  medical  procedures.  In  this  paper,  large 
deformation  is  considered  in  the  FE  formulation,  and  soft 
tissues  are  modeled  as  incompressible  hyperelastic  materials 
based  on  the  experimental  data  of  [9],  Viscoelastic  response 
is  not  taken  into  account  in  this  model  as  we  do  not  expect 
it  to  be  important  on  the  short  time  scales,  such  as  in  breast 
medical  applications.  Nevertheless,  it  is  possible  to  incorporate 
viscoelasticity  into  our  model  provided  that  relevant  parameters 
are  available.  Boundary  conditions  in  medical  applications 
are  usually  displacement  type  which  are  known  relatively 
accurately.  In  multimodality  imaging  which  involves  X-ray 
mammography,  the  boundary  conditions  are  not  fully  known 
as  the  contact  surface  cannot  be  identified  easily.  However,  a 
contact  problem  formulation  as  presented  in  [1]  can  be  used 
which  only  requires  the  magnitude  of  the  compression  plate 
movement. 

For  breast  FE  mesh  generation,  marching-cubes-based  tech¬ 
niques  such  as  [10]  can  be  implemented.  However,  the  tetra¬ 
hedral  elements  used  in  these  techniques  are  well  known  to 
have  unfavorable  characteristics,  such  as  slow  convergence  with 
mesh  refinement.  Also,  these  elements  may  exhibit  overstiff¬ 
ening  and  volumetric  locking  especially  with  incompressible 
material  [11],  [12].  Accordingly,  with  biological  soft  tissues 
which  are  known  to  be  incompressible,  using  hexahedral  ele¬ 
ments  in  the  FE  model  is  cost  effective  and  more  accurate.  To 
our  knowledge,  there  has  been  no  work  done  for  patient-specific 
mesh  generation  of  the  breast  using  hexahedral  elements.  How¬ 
ever,  computed  tomography  (CT)  has  been  used  by  a  number 
of  investigators  to  automate  mesh  generation  in  bone.  Many 
investigators  [2],  [13]-[17]  have  applied  edge-detection  algo¬ 
rithms  on  CT  images  to  generate  boundary  contours  of  the  struc- 
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tures  in  each  image,  then  employed  some  techniques  to  gen¬ 
erate  the  internal  mesh.  This  strategy  often  requires  a  signifi¬ 
cant  investment  of  time  by  an  expert  user,  and  is  incapable  of 
automatically  generating  reasonably  accurate  meshes  for  het¬ 
erogeneous  and  geometrically  complex  objects.  Other  investiga¬ 
tors  [  1 8]— [20]  developed  techniques  that  directly  convert  image 
voxels  to  eight  noded  hexahedral  elements.  The  resulting  ob¬ 
ject  surfaces  in  these  techniques  arc  characterized  by  abrupt 
transitions  and  right  angles.  As  a  result,  FEM  leads  to  less  ac¬ 
curate  results,  especially  on  the  surface  [21],  [22].  This  paper 
presents  two  methods  for  generating  three-dimensional  (3-D) 
meshes  of  the  breast  using  magnetic  resonance  imaging  (MRI) 
data.  The  first  method  converts  image  voxels  to  hexahedral  ele¬ 
ments  and  is  similar  to  the  said  methods  except  that  edge  detec¬ 
tion  is  done  semi-automatically  via  segmentation.  The  second 
method  is  based  on  the  transfinite  interpolation  (TI)  technique 
[23]  which  is  more  accurate  and  efficient  in  creating  smooth  sur¬ 
faces. 

To  validate  the  mesh  generation  component  of  this  model, 
a  cubic  agarose  phantom  with  a  cylindrical  inclusion  was 
constructed  and  imaged  using  MRI.  Using  these  images,  FE 
meshes  were  created  using  the  meshing  techniques  presented 
in  this  paper.  To  assess  the  accuracy  of  these  methods,  another 
mesh  was  created  manually  based  on  the  known  geometry 
of  the  phantom.  The  FEM  was  applied  to  these  meshes  to 
calculate  the  deformation  and  stress  distribution  resulting 
from  a  lateral  compression  similar  to  X-ray  mammography 
procedures.  The  results  obtained  from  both  meshes  compare 
well  with  the  manually  created  mesh  results  except  for  the 
inclusion’s  circular  interface  where,  as  expected,  the  Tl-based 
mesh  produced  better  results.  To  demonstrate  a  potential 
clinical  application  of  this  model,  MR  images  of  a  healthy 
volunteer’s  breast  was  acquired  and  a  FE  mesh  was  created 
using  the  meshing  techniques.  A  lateral  compression  similar 
to  X-ray  mammography  procedures  was  simulated  using  the 
Tl-based  mesh,  and  calculated  MR  images  are  presented  in  the 
postcompression  configuration. 


II.  Methods 


A.  Governing  Equations 

We  use  a  3-D  FE  model  to  predict  breast  tissue  deformation 
based  on  the  biomechanical  parameters  of  the  breast  tissues. 
While  undergoing  medical  procedures,  the  breast  often  deforms 
significantly.  Accordingly,  linear  elasticity  with  the  infinites¬ 
imal  deformation  formulation  is  not  appropriate  to  formulate 
the  FE  model.  As  a  result,  we  used  a  finite  deformation  for¬ 
mulation  in  conjunction  with  hyperelastic  material.  In  a  static 
condition,  the  equilibrium  equation  of  a  body  under  externally 
applied  forces  is 


dan  ,  dai2  ,  0ai3  ,  t  _  n 
»  rk  i  o  i  Ji  — 


*  =  1,2,3  (1) 


dx\  dx2  dxs 
where  the  as  denote  the  Cauchy  (true)  stress  tensor  components 
and  fi  represents  the  body  force.  In  the  finite  deformation  for¬ 
mulation,  the  current  position  of  a  point  x  is  obtained  by  adding 
the  displacement  u  to  the  corresponding  reference  position  of 
the  same  point  X,  i.e., 

x  =  X  4-  u.  (2) 


The  stresses  are  related  to  displacements  (u)  via  a  strain  energy 
function  U  which  is  defined  as  a  function  of  the  strain  invariants. 
These  invariants  arc  obtained  based  on  the  deformation  gradient 
F  as  follows: 


F  = 


Ox 

OX' 


(3) 


To  separate  the  deviatoric  and  volumetric  effects,  F  =  J-1/3F 
is  defined,  where  J  =  det(F)  is  the  third  strain  invariant  which 
represents  the  total  volume  change  at  the  point.  Using  B  = 
F.F  ,  the  other  strain  invariants  /1?  /2  can  be  defined  as  fol¬ 
lows: 


h  =tr(B) 

72  =  i(7?-tr(B.B)). 

Assuming  that  a  strain  energy  function  U  =  U(I i,  I 2,  J)  is  se¬ 
lected,  the  true  stress  tensor  can  be  obtained  using  the  following 
equation  [11]: 


a  =  i 
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where  I  is  the  unit  matrix.  Considering  the  boundary  conditions, 
(1)  and  (4)  must  be  solved  simultaneously  to  calculate  the  dis¬ 
placements  u.  These  equations  are  solved  numerically  through 
the  FE  method.  In  this  method,  after  the  domain  is  discretized 
into  a  number  of  homogeneous  elements,  a  nonlinear  system 
of  equations  is  formed  that  is  solved  iteratively  using  Newton’s 
method.  This  formulation  is  implemented  in  ABAQUS,  a  com¬ 
mercial  FE  software  [11].  Discretization,  i.e.,  mesh  generation, 
especially  in  medical  applications,  is  the  most  tedious  step  in  a 
FE  analysis.  This  step  is  described  in  the  following  sections. 


B .  Image  Acquisition 

We  have  developed  meshing  techniques  which  use  MR 
images  as  an  input  to  create  FE  meshes.  In  this  work,  two  sets 
of  MR  images  were  acquired  for  a  phantom  and  a  volunteer’s 
breast.  The  MR  images  were  acquired  using  a  GE  SIGNA 
1.5-T  scanner  (GE  Medical  Systems,  Milwakee,  WI).  For  the 
phantom,  a  5.0-in  surface  coil  was  used  to  acquire  Tl -weighted 
axial  images.  A  two-dimensional  (2-D)  spin-echo  pulse  se¬ 
quence  was  used  with  90°  flip  angle,  TR/TE  =  400.0/10.5 
msec  timing  properties,  12x12  cm  field  of  view  (FOV),  256  x 
128  resolution  and  1 .0-mm  slice  thickness.  A  magnitude  image 
of  this  phantom  is  shown  in  Fig.  1 ,  in  which  the  inclusion  can  be 
seen  with  subtle  contrast.  For  the  breast,  Tl  weighted  sagittal 
images  were  obtained  using  body  coil.  A  2-D  spin-echo  pulse 
sequence  was  used  with  90°  flip  angle,  TR/TE  =  300.0/9.0 
ms  timing  properties,  20  x  20  cm  FOV,  256  x  128  resolution, 
3.0-mm  slice  thickness,  and  three  acquisitions. 

C.  Image  Segmentation  and  Edge  Detection 

The  purpose  of  segmentation  is  to  separate  different  mate¬ 
rials  within  the  object  of  interest.  This  is  done  based  on  the  MR 
signal  intensity  of  different  tissues.  In  the  breast,  the  two  main 
components,  i.e.,  adipose  and  fibroglandular  tissues,  must  be 
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(a)  (b) 

Fig.  1 .  (a)  A  typical  MR  image  and  (b)  a  corresponding  segmented  image  of  the  phantom. 


segmented.  For  edge  detection,  as  the  breast  tissues  have  sig¬ 
nificant  signal  intensity  compared  to  the  background  with  only 
noise  signal,  segmentation  results  in  automatic  breast  edge  de¬ 
tection.  For  segmentation,  3-D  MRI  data  of  the  breast  is  con¬ 
verted  into  AnalyzeAVW  2.5  [24]  database  where  a  standard 
thresholding-based  tool  was  used  to  yield  segmented  images 
that  are  suitable  for  FE  mesh  generation. 

D.  Voxel-Based  (VB)  Mesh  Generation  Technique 

Because  of  the  complex  geometry  of  the  breast  and  the 
structure  of  the  fibroglandular  tissue,  common  mesh  generators 
require  a  significant  amount  of  user  interaction  which  is  very 
time  consuming.  As  an  attempt  to  minimize  user  interaction 
throughout  the  mesh  generation  step,  we  have  developed 
an  automatic  VB  mesh  generator  with  the  aim  of  creating 
relatively  accurate  FE  meshes  of  complex  volumetric  data 
efficiently.  In  fact,  the  only  step  which  requires  user  interaction 
in  this  technique  is  the  segmentation  step.  To  construct  the 
mesh,  eight  noded  hexahedral  elements  are  chosen.  This  type  of 
element  is  well  known  to  have  many  favorable  characteristics 
over  tetrahedral  elements  which  can  demonstrate  overstiffening 
and  locking  behavior  [12].  Right-angled  hexahedrals,  however, 
cannot  model  curved  surface  geometries  which  in  turn  leads  to 
some  inaccuracies  in  surface  strains.  To  improve  the  mesh  in 
this  respect,  we  used  the  smoothing  technique  of  Camacho  et 
al.  [25]  which  smooths  irregular  boundaries  at  model  surfaces 
and  material  interfaces. 

In  the  VB  methods,  since  one  FE  is  generated  corresponding 
to  each  voxel,  high  resolution  MR  images  may  lead  to  huge 
FE  meshes  which,  due  to  current  computer  power  limitations, 
cannot  be  processed.  Moreover,  beyond  a  certain  point,  finer  FE 
meshes  do  not  practically  improve  the  solution  accuracy.  There¬ 
fore,  reducing  image  resolution  is  required  to  meet  computa¬ 
tional  pbwer  constraints.  For  this  purpose,  a  larger  voxel  size 
is  selected  and  using  the  segmented  images,  the  adipose  and  fi¬ 
broglandular  tissues  voxels  within  each  new  voxel  are  counted. 


The  new  voxel  will  then  be  assigned  to  be  as  the  one  with  the 
larger  count. 

After  segmentation  and  resolution  reduction,  the  mesh  gener¬ 
ation  process  is  followed  by  node,  and  later  element  generation. 
For  node  generation,  based  on  the  pixel  size,  slice  thickness  and 
the  coordinate  system  of  the  MR  scanner,  the  nodes  are  arranged 
in  the  form  of  a  3-D  rectangular  lattice  which  confines  the  breast 
volume.  The  node  coordinates  of  the  lattice  are  calculated  and 
placed  in  an  array.  Working  with  the  regular  lattice,  later,  facil¬ 
itates  element  generation  significantly  because  each  neighbor 
node  can  be  determined  easily  and  accessed  very  quickly  in  the 
node  array.  The  lattice,  however,  represents  regions  both  inside 
and  outside  the  breast.  Regions  outside  the  breast  are  determined 
during  the  element  generation  process  and  then  discarded  from 
the  nodes  array  after  element  generation  is  done. 

For  element  generation,  the  segmented  image  of  each  slice 
is  loaded  by  the  computer  and  stored  in  a  2-D  array.  To  facil¬ 
itate  finding  the  boundary  of  each  slice  automatically,  we  as¬ 
sume  that  the  fibroglandular  tissue  is  surrounded  by  the  adi¬ 
pose  tissue  portion.  If  this  is  not  the  case,  a  very  thin  one  pixel 
size  layer  which  represents  adipose  tissue  can  be  added  on  the 
image  manually  using  XV  3.10a  image  processor  [26]  to  close 
the  boundary  where  necessary.  The  boundary  points  are  found 
by  searching  the  image  array  line  by  line  first  from  left  to  right 
and  then  from  right  to  left  until  an  adipose  tissue  pixel  on  the 
left  and  right  boundaries  is  reached.  These  points  are  stored  in 
two  arrays  that  are  later  used  in  element  generation.  Elements 
corresponding  to  adipose  and  fibroglandular  tissues  are  gener¬ 
ated  in  each  slice  row  by  row  starting  from  the  left  boundary 
point  as  the  local  first  element  node.  Other  nodes  in  each  ele¬ 
ment  are  determined  based  on  the  local  numbering  scheme  of 
hexahedral  element.  Element  generation  in  a  row  is  continued 
until  the  point  prior  to  the  right  boundary  point  is  reached.  This 
process  is  repeated  for  each  row  to  obtain  a  slice  mesh  and  then 
for  each  slice  to  cover  the  entire  breast.  Slice  meshes  are  then 
stacked  to  reconstruct  a  complete  3-D  mesh.  Skin  elements  are 
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Fig.  2.  (a)  A  unit  square  (logical  space),  (b)  A  typical  physical  space,  (c)  Breast  boundaries  in  the  TI  meshing  technique. 


generated  using  four  noded  quadrilateral  membrane  elements. 
These  elements  are  generated  by  identifying  the  outer  face  of 
elements  on  the  breast  surface.  Finally,  for  smoothing,  the  node 
coordinates  are  modified  using  the  smoothing  technique  of  Ca¬ 
macho  et  ai  [25].  In  this  technique,  a  smoothing  coefficient 
which  varies  from  zero  for  no  smoothing  to  0.5  for  maximum 
smoothing,  must  be  selected  by  the  user.  While  higher  factors 
result  in  smoother  surfaces  and  interfaces,  i.e.,  better  geometry 
modeling,  element  distortion  caused  by  this  process  leads  to  less 
accurate  FE  analysis.  Therefore,  to  achieve  optimal  results,  the 
smoothing  coefficient  must  be  selected  or  determined  carefully. 
The  meshing  algorithm  is  implemented  by  creating  a  MATLAB 
program  which  inputs  the  segmented  axial,  coronal  or  sagittal 
images  obtained  from  MRI  data,  generates  the  FE  mesh,  and 
outputs  an  input  file  compatible  with  ABAQUS  FE  software 
[11].  Henceforth,  this  meshing  program  will  be  referred  to  as 
VBMESH. 

E.  Transfinite  Interpolation  Mesh  Generation  Technique 

This  technique  is  based  on  the  idea  of  mapping  a  unit  square 
(logical  space)  to  any  given  shape  enclosed  by  four  boundaries 
(physical  space)  [23].  To  implement  this  idea  in  breast  FE 
meshing,  after  image  segmentation,  the  breast  boundary  in 
each  slice  is  divided  into  three  segments.  Adding  the  boundary 
representing  the  chest  wall,  as  shown  in  Fig.  2,  the  breast  area 
will  be  enclosed  by  four  boundaries.  Although  not  necessary, 
these  boundaries  can  be  fitted  by  polynomials  to  smooth  out  the 
surface  which  is  corrupted  as  a  result  of  segmentation  errors. 
After  obtaining  the  breast  boundaries  in  each  slice  image,  the 
nodes  are  generated  by  mapping  the  nodes  of  the  unit  square 
grid  to  a  set  of  nodes  distributed  inside  the  boundaries  using 
the  following  equation: 

X(t,  7?)  =  (1  -  V)xb(0  +  vxt(0  +  (1  -  ox,(v) 

+  £Xr{v)  -  SvXti i)  -  £(i  -  v)Xt(i) 

-  „(1  -  0*t(0)  -  (1  -  0(1  -  v)xb(0)  (5) 

where  £  and  rj  represent  the  unit  square  variables  and  Xt,  Xt, 
Xi  and  Xr  represent  the  breast’s  bottom,  top,  left  and  right 
boundaries  respectively.  Element  connections  determination  in 
this  technique  is  straight  forward  as  it  is  done  using  the  nodes 
of  the  unit  square’s  grid.  Using  the  same  logic,  the  FE  mesh  of 
the  skin  is  created  by  finding  node  connections  of  three  of  the 


unit  square’s  circumferential  lines  corresponding  to  the  breast’s 
surface.  After  elements  in  each  slice  are  generated,  as  in  the  pre¬ 
vious  method,  they  are  stacked  to  construct  a  3-D  mesh. 

Compared  to  the  VB  technique,  this  method  requires  a  signif¬ 
icantly  smaller  number  of  elements  to  represent  the  true  geom¬ 
etry  of  the  breast’s  surface.  Furthermore,  this  method  is  capable 
of  locally  refining  the  mesh  at  an  ROI  by  choosing  smaller  in¬ 
tervals  in  the  corresponding  logical  space.  Similar  to  VBMESH, 
this  algorithm  was  implemented  by  a  MATLAB  program  which 
will  be,  henceforth,  referred  to  as  TIMESH. 

III.  Results 

A.  Meshing  Techniques  Validation  and  Comparison 

To  validate  the  meshing  techniques,  the  agarose  phantom  il¬ 
lustrated  in  Fig.  3  was  constructed.  This  phantom  consists  of  a 
block  of  agar  (2%)  with  a  hard  cylindrical  inclusion  made  of  a 
mixture  of  2%  agar  and  1 0%  glass  beads  centered  at  (32.8, 30.5). 
The  glass  beads  were  used  to  enhance  the  contrast  in  the  MR 
image  for  segmentation  purposes.  We  assume  that  the  Young’s 
moduli  of  the  block  and  the  inclusion  are  respectively  1 1 .0  Kpa 
and  55.0  Kpa.  The  block  is  also  assumed  to  undergo  a  static 
loading  due  to  9.0-mm  compression  resulting  from  a  laterally 
moving  rigid  plate  toward  another  stationary  plate  as  illustrated 
in  Fig.  3.  The  phantom  was  scanned  using  a  5.0-in  surface  coil 
and  a  set  of  60  MR  images  was  acquired.  These  images  were 
segmented  and  the  resolution  was  reduced  by  a  factor  of  four  in 
the  x  and  y  directions  and  to  16  slices  in  the  z  direction.  An  MR 
image  of  the  phantom  and  a  corresponding  segmented  image 
are  shown  in  Fig.  1.  The  obtained  images  were  processed  by 
VBMESH  using  a  smoothing  coefficient  of  0.3  and  a  FE  mesh 
was  created.  By  dividing  the  inclusion’s  circumference  into  four 
segments  and  fitting  each  segment  with  a  quadratic  polynomial, 
another  mesh  was  created  using  TIMESH.  Finally,  based  on  the 
known  geometry  of  the  phantom,  a  mesh  was  created  manu¬ 
ally  as  a  bench  mark.  Corresponding  to  each  mesh,  an  input  file 
compatible  with  ABAQUS  was  created.  These  input  files  were 
preprocessed  and  their  corresponding  meshes  are  displayed  in 
Fig.  4.  Based  on  the  nonlinear  finite  deformation  FE  formula¬ 
tion,  the  meshes  were  analyzed  using  ABAQUS  5.8.1 ,  and  as  a 
result,  the  displacements,  strains  and  stresses  were  calculated.  It 
is  well  known  that  in  FE  analysis  the  most  significant  errors  are 
encountered  at  material  interfaces  and  curved  surfaces  which 
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66.0  mm 


Fig.  3.  A  schematic  of  the  phantom  and  its  compression  model. 


(a)  (b)  (c) 

Fig.  4.  (a)  VB  FE  mesh,  (b)  Tl-based  FE  mesh,  and  (c)  A  manually  created  FE  mesh  of  the  phantom. 


are  not  represented  appropriately  by  the  FE  mesh  [22].  There¬ 
fore,  to  assess  the  accuracy  of  the  presented  meshing  techniques, 
we  will  focus  on  the  nodal  values  at  the  inclusion’s  circumfer¬ 
ence.  As  an  example,  the  results  of  a  displacement,  a  strain,  and 
a  stress  component  at  this  region  in  the  middle  plane  are  dis¬ 
played  in  Fig.  5,  which  qualitatively  indicates  the  better  perfor¬ 
mance  of  the  Tl-based  technique.  To  make  a  quantitative  com¬ 
parison,  the  average  E  and  standard  deviation  a  of  the  relative 
error  of  the  displacement,  strain  and  stress  over  the  inclusion’s 
circumferential  nodes  are  calculated  and  summarized  in  Table  I. 
The  results  clearly  indicate  that  the  Tl-based  technique  leads 
to  more  accurate  FE  meshes.  Moreover,  compared  to  the  dis¬ 
placements  errors,  the  strain  and  stress  errors  are  significantly 
larger.  This  can  be  justified  based  on  the  fact  that  calculating 


both  strains  and  stresses  involves  calculating  the  derivatives  of 
displacements  which  leads  to  error  magnification. 

B.  Breast  Image  Nonrigid  Registration 

To  examine  our  model  in  a  complex  clinical  problem, 
we  applied  it  to  a  breast  image  registration  problem  where 
the  breast  was  compressed  by  two  plates  causing  nonrigid 
deformation.  For  this  purpose,  sagittal  MR  images  of  a  breast 
of  a  healthy  volunteer  were  acquired  at  6-mm  intervals  at  a 
resolution  of  0.625  mm/pixel.  These  images  were  segmented 
with  Analyze AVW  2.5  and  the  resulting  images  were  processed 
to  create  FE  meshes.  A  typical  sagittal  MR  image  of  the  breast 
and  a  corresponding  segmented  image  are  shown  in  Fig.  6(a) 
and  (b).  The  objective  here  was  to  calculate  MR  images  of 
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TABLE  1 

Avkragk  Rhi.ativk  Errors  E  and  Standard  Dhviation  ct  or 
Displachmhnts,  Strains,  and  Strkssks  okthk  Cai.cui.athi>  Mhshhs  at 

THK  INCLUSION’S  SURFACH 


Displacements 

U 1 

U2 

u3 

- 

- 

- 

VB  Mesh:  E  (%) 

7.50 

7.93 

5.67 

- 

- 

- 

<7  (%) 

1.38 

4.07 

3.07 

- 

- 

- 

TI  Mesh:  E  (%) 

1.70 

2.04 

0.84 

- 

- 

- 

<7(%) 

0.50 

1.69 

0.78 

- 

- 

- 

Strains 

Cll 

^22 

£33 

£12 

Cl  3 

^23 

VB  Mesh:  E  (%) 

18.14 

16.94 

6.32 

16.56 

20.86 

20.89 

<7  (%) 

7.00 

9.79 

4.79 

11.00 

15.23 

20.91 

TI  Mesh:  E  (%) 

14.24 

8.08 

9.82 

17.89 

13.10 

16.33 

<7(%) 

11.25 

7.17 

2.83 

12.83 

10.21 

17.27 

Stresses 

<711 

CT22 

<7.33 

<712 

<713 

<723 

VB  Mesh:  E  (%) 

17.15 

14.55 

9.37 

16.56 

20.86 

20.89 

a  (%) 

9.32 

14.89 

10.32 

11.00 

15.23 

20.91 

TI  Mesh:  E  (%) 

13.20 

8.21 

12.03 

17.89 

13.10 

16.33 

a  (%) 

10.70 

8.31 

5.36 

12.83 

10.21 

17.27 

this  breast  as  if  it  underwent  8.0~mm  compression.  For  this 
purpose,  two  FE  meshes  were  first  created  using  the  VB  and 
Tl-based  techniques.  After  image  segmentation  and  resolution 
reduction  by  a  factor  of  4  in  each  direction,  the  resulting  images 
were  processed  by  VBMESH  using  a  smoothing  coefficient 
of  0.3.  The  resulting  mesh  consists  of  16841  elements  and 
15  939  nodes  and  is  characterized  by  an  abrupt  surface.  Using 
the  segmented  images,  the  Tl-based  mesh  was  created  using 
TIMESH.  This  mesh  was  calculated  by  fitting  a  third-order 
polynomial  to  the  boundary  representing  the  chest  wall  while 
fitting  quadratic  polynomials  to  the  breast  surface’s  three 
segments.  The  resulting  mesh  consists  of  2280  elements  and 
2059  nodes  and,  unlike  the  VB  mesh,  has  a  smooth  surface.  The 
FE  meshes  of  one  slice  created  using  the  VB  and  the  Tl-based 
techniques  are  depicted  in  Fig.  6(c)  and  (d).  Fig.  7  depicts 
the  VB  mesh,  the  Tl-based  mesh  of  the  entire  breast,  and  the 
breast  FE  model  undergoing  compression.  The  adipose  and 
fibroglandular  tissues  were  assumed  to  be  incompressible  and 
hyperelastic,  and  the  corresponding  parameters  were  obtained 
by  fitting  experimental  data  [9]  to  a  hyperelastic  Neo-Hookean 
model.  For  this  purpose,  the  data  represented  as  Young’s 


modulus  versus  strain  were  converted  to  stress  versus  strain  by 
first  fitting  quadratic  and  third-order  polynomials  to  the  adipose 
and  fibroglandular  tissues  data,  respectively,  then  integrating 
over  the  strain.  The  resulting  fitting  polynomials  of  the  Young’s 
modulus  versus  strain  of  the  adipose  and  fibroglandular  tissues 
are  as  follows: 

E  =  0.5197e2  +  0.0024r  +  0.0049;  for  fat 

E  =  123.8889r3  -  11.7667e2 

+().G9G9c  +  0.0121;  for  parenchyma. 

To  the  authors’  knowledge,  the  experimental  data  [9]  arc  the 
only  available  breast  tissue  hyperelastic  properties  data  obtained 
based  on  the  large  deformation  theory  considerations.  Neverthe¬ 
less,  as  the  standard  deviation  of  this  data  is  high,  more  accurate 
and  reliable  data  is  required  to  produce  better  image  registra¬ 
tion  results.  Skin,  in  general,  is  well  known  to  be  anisotropic 
and  hyperelastic  [27].  However,  it  can  be  approximated  by  an 
isotropic  and  linear  elastic  model  provided  that  the  strain  does 
not  exceed  50%  [27],  [28].  Accordingly,  as  the  skin  strain  in 
this  application  is  expected  to  be  well  below  50%,  it  was  as¬ 
sumed  to  be  elastic  with  a  Young’s  modulus  of  10.0  kPa  [28] 
and  a  thickness  of  1.0  mm.  The  chest  wall  was  assumed  to  pro¬ 
vide  zero  displacement  boundary  conditions  while  the  specified 
displacement  or  force  boundary  conditions  on  the  breast’s  sur¬ 
face  is  not  known  because  the  contact  surface  between  the  breast 
and  the  plates  increases  over  the  course  of  compression.  Accord¬ 
ingly,  this  problem  was  formulated  as  a  3-D  contact  problem  and 
ABAQUS’  contact  solver  with  the  finite  deformation  theory  for¬ 
mulation  was  used.  This  highly  nonlinear  problem  was  solved 
iteratively  and,  as  a  result,  the  breast  tissue  displacements  were 
calculated.  These  displacements,  in  combination  with  the  MR 
images  of  the  uncompressed  breast,  were  used  to  obtain  simu¬ 
lated  MR  images  of  the  compressed  breast.  A  number  of  con¬ 
secutive  sagittal  MR  images  of  the  breast  before  compression, 
along  with  their  corresponding  simulated  images  after  compres¬ 
sion,  are  depicted  in  Fig.  8.  In  this  figure,  reasonable  translation 
and  tissue  expansion  can  be  observed  by  comparing  the  two  sets 
of  images.  For  example,  image  (g')  from  the  simulated  set  and 
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Fig.  6.  (a)  A  MR  image  and  (b)  a  corresponding  segmented  image,  (c)  VB  FE  mesh,  and  (d)  Tl-based  FE  mesh  of  one  breast  slice. 


image  (f)  from  the  MR  image  set,  which  approximately  cor¬ 
respond  to  the  same  slice,  depict  similar  anatomical  features. 
The  same  can  be  observed  in  other  image  pairs  such  as  (a,  6'), 
(b:  c'),  etc. 

IV.  Discussion 

We  have  presented  a  FE  model,  based  on  biomechanical 
principles,  to  predict  breast  tissue  deformations.  This  model 
can* serve  as  a  guideline  in  numerous  clinical  applications,  such 
as  breast  image  registration,  multimodality  data  fusion,  breast 
surgery,  and  biopsy.  The  first  component  of  this  model  is  FE 


meshing,  for  which  two  meshing  techniques  were  presented 
and  evaluated.  These  techniques,  which  require  minimal 
user  interaction,  use  MRI  images  of  the  breast  to  produce 
patient- specific  FE  meshes.  Using  a  numerical  experiment,  the 
Tl-based  meshing  technique,  which  requires  fewer  elements  to 
represent  the  geometry  of  curved  surfaces  relatively  accurately, 
proved  to  generate  more  accurate  displacements  and  stresses. 
Moreover,  unlike  the  VB  technique,  this  technique  is  capable  of 
refining  the  mesh  at  an  arbitrary  ROI.  The  adipose  and  fibrog- 
landular  tissues  are  modeled  as  incompressible  and  hyperelastic 
materials  while  the  skin  is  modeled  using  membrane  elements 
with  linear  elastic  properties.  These  tissues  are  assumed  to 
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undergo  large  deformations.  As  for  boundary  conditions,  the 
'chest  wall  region  is  assumed  to  provide  zero  displacements 
boundary  conditions  while  the  specified  boundary  conditions 
depend  on  the  clinical  application.  For  example,  in  the  breast 
image  registration  problem,  as  was  previously  presented, 
this  boundary  condition  type  is  not  known  explicitly  and  the 
problem  is  formulated  as  a  3-D  contact  problem.  An  image 
registration  example  was  presented  to  evaluate  the  FE  model’s 
effectiveness  in  complex  clinical  applications.  The  results 
summarized  in  Fig.  8  qualitatively  prove  the  merits  of  this 
model.  Quantitative  validation,  however,  is  beyond  the  scope  of 
this  paper  as  it  requires  more  reliable  breast  tissue  hyperelastic 
properties  data.  Research  is  underway  in  our  laboratory  to 
acquire  such  data  which  will  pave  the  way  for  conducting  this 
validation. 
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